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Abstract 

A new damage model based on a micromechanical analysis of cracked [±0/9O n ] s 
laminates subjected to multiaxial loads is proposed. The model predicts the onset 
and accumulation of transverse matrix cracks in uniformly stressed laminates, the 
effect of matrix cracks on the stiffness of the laminate, as well as the ultimate failure 
of the laminate. The model also accounts for the effect of the ply thickness on the 
ply strength. Predictions relating the elastic properties of several laminates and 
multiaxial loads are presented. 


Key words: Failure criteria, Micromechanics, Damage Mechanics. 


1 Introduction 


The aerospace industry is committed to improve the performance of aircraft 
whilst reducing emissions and weight. Such a goal can be achieved by the use 
of advanced polymer-based composite materials, that have excellent proper- 
ties for aerospace applications, such as low density, and fatigue and corrosion 
resistance. 

The design procedure used for advanced composite structures relies on a 
’building-block’ approach [1], where a large number of experimental tests are 
performed throughout the product development process. The use of improved 
analytical or numerical models in the prediction of the mechanical behavior 
of composite structures can significantly reduce the cost of such structures. 
Such models should predict the onset of material degradation, the effect of the 



non-critical damage mechanisms on the stiffness of the laminate, and ultimate 
structural failure. 

Strength-based failure criteria are commonly used to predict failure in com- 
posite materials [2]-[7]. Failure criteria predict the onset of the several damage 
mechanisms occurring in composites and, depending on the laminate, geome- 
try and loading conditions, may also predict structural collapse. 

In multidirectional composite laminates, damage accumulates during the load- 
ing process. Final failure occurs as a result of damage accumulation and stress 
re-distribution. The ultimate failure load is higher than the initial damage 
predicted by strength-based failure criteria. Furthermore, stress- or strain- 
based failure criteria cannot represent size effects that occur in quasi-brittle 
materials [8]. 

Simplified models, such as the ply discount method where some scalar compo- 
nents of the stiffness tensor are reduced to approximately zero when damage is 
predicted, are often used by the industry to predict ultimate laminate failure. 
However, these methods cannot represent with satisfactory accuracy the pro- 
gressive reduction of the stiffness of a laminate as a result of the accumulation 
of matrix cracks. 

Constitutive laws based on Continuum Damage Mechanics (CDM) have been 
proposed to predict the material response, from the onset of damage up to final 
collapse [9]- [14]. Although the existing CDM models can accurately predict the 
evolution of damage, they usually rely on fitting parameters that need to be 
measured at laminate level, such as the critical values of thermodynamic forces 

[14]- 

Three-dimensional models based on CDM use material properties measured 
at ply level. Damage mechanisms such as transverse matrix cracks are rep- 
resented by strain-softening constitutive models and bands of localized de- 
formation. The implementation of strain-softening models in finite element 
models may cause convergence difficulties during the iterative solution proce- 
dure. Furthermore, strain-softening models require regularization techniques 
to provide mesh-independent solutions [15,16]. 

Alternative methods based on the combination of clastic analysis of cracked 
plies and finite Fracture Mechanics provide the basis for an accurate repre- 
sentation of the response of composite materials [17]- [30]. Micromechanical 
models have been developed to predict the initiation and evolution of trans- 
verse matrix cracks under either in-plane shear or transverse tension. Gener- 
alizations of these models are required for the more usual case of multiaxial 
loading. 
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In order to predict ultimate failure, a micromechanical model needs to be used 
in combination with a fiber failure criterion. Furthermore, a general method- 
ology to predict laminate strength should be able to predict matrix cracking 
under high values of transverse compression. These damage mechanisms usu- 
ally correspond to the final failure of laminates under uniform states of stress 
[ 6 ], 

The objective of this work is to define a new damage model based on microme- 
chanical models of transverse matrix cracks to predict the onset and evolution 
of transverse matrix cracks under multiaxial loading. A new constitutive model 
is derived based on the thermodynamics of irreversible processes. Continuum 
Damage Mechanics provides a rigorous framework to define the constitutive 
model and to develop the corresponding computational implementation. The 
model proposed herein represents transverse matrix cracks as distributed dam- 
age and does not require any regularization to yield mesh-independent results 
when implemented in finite element models. Furthermore, the model is able 
to predict the onset and propagation of matrix transverse cracks under multi- 
axial loading as well as the final failure of uniformly stressed laminates where 
a periodic distribution of transverse matrix cracks can be assumed. 


2 Micromechanical model of transverse matrix cracks under mul- 
tiaxial loading 


The proposed continuum damage model is based on two major components: a 
set of stress based failure criteria and a micromechanical model of transverse 
matrix cracking in multidirectional laminates. 

The failure criteria define the onset of transverse matrix cracking, i.e. the 
activation of the damage variables. A micromechanical model of transverse 
matrix cracks is required to define the evolution of the damage variables. 

Several micromechanical models of transverse matrix cracks that have been 
proposed in the literature [19,30] can be used within the framework devel- 
oped here. The micromechanical model proposed by Tan and Nuismer [17,18] 
accounts for the effects of the adjoining plies on the homogenized elastic prop- 
erties of a cracked ply. This micromechanical model is the basis of the develop- 
ments presented in this work. The models are considered as micromechanical 
in the sense that they are based on a boundary value problem that explicitly 
represents transverse matrix cracks. However, the fibers and the matrix are 
not represented explicitly. 
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Using the assumption of generalized plane strain, Tan and Nuismer [17,18] 
developed a model able to relate the density of transverse matrix cracks in a 
central 90° ply to the homogenized clastic properties of that ply. The model 
developed by Tan and Nuismer was used for the prediction of the evolution 
of transverse matrix cracks under either in-plane shear or transverse tensile 
stresses. 

The laminates under investigation are symmetric and balanced with a [±#/90 n ] s 
layup containing a periodic distribution of transverse matrix cracks, as shown 
in Figure 1. The micromechanical analysis of a balanced symmetrical laminate 
requires the division of the laminate in two sub-laminates: the 90° layers in 
the middle layer (sublaminate 1), and the outer plies (sublaminate 2). 



Fig. 1. One-quarter of the laminate. 

Sublaminate 1 is taken as transversely isotropic. The outer layers, sublami- 
nate 2, may contain several layers with different fiber orientations but must 
be orthotropic in the laminate global system ( x , y and z coordinates). The 
expressions for the compliance and stiffness tensors for the laminate, sublam- 
inate 1, and sublaminate 2 are shown in Appendix A. In what follows, the 
superscripts (1) and (2) corresponds to the sublaminate 1 and 2, respectively. 
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2.1 Damage dependent constitutive tensor 


The stiffness tensor A of the balanced and symmetric laminate shown in Figure 
1 can be written as a function of one-half of the distance between transverse 
matrix cracks (L) in sublaminate 1: 


ML) 


A(L ) n A(L) l2 0 

A{L) 21 A(L) 22 0 

o o A(L) m 


(i) 


where the terms Aij(L),i,j = 1,2,6 are obtained from the Tan and Nuismer 
model [18] which is summarized in Appendix A. 


The undamaged stiffness matrix Q of the laminate is obtained from lamination 
theory using the undamaged stiffness matrix of sublaminate 1, QW, and the 
stiffness matrix of sublaminate 2, Q( 2 ). Therefore, 


( 5 < 5 « - h(2> Q ( ff) ( 2 ) 

where h = h w + hM . 

Assuming that the degradation due to the transverse matrix cracks only occur 
in sublaminate 1, the damaged stiffness tensor of laminate 1 is given as: 


4f(i) = ^5 (hA v (L) - h^Q\f) (3) 

where A^(L) ( i,j = 1,2,6) are the scalar components of the stiffness tensor 
of sublaminate 1. These components are a function of the distance between 
transverse matrix cracks (L, defined in Figure 1). 

Having defined A^(L), it is possible to calculate the effective transverse mod- 
ulus, Poisson ratio, and shear modulus of the 90° ply: 


E?(L) = Eg\L) = iff(L) 


^22 U) 



a { 3(l) 

4 * 2 \l) 


g ( 2(l) = &»(l) = 4«'(i) 


( 4 ) 

( 5 ) 

( 6 ) 
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One noticeable feature of this model is that the degradation of the transverse 
modulus is equal to the degradation of the minor Poisson ratio: 


V 2 l(L) = V21 , . 

e£\l) e 2 

Therefore, the quotient V 21 /E 2 is not a function of damage. This observation 
is in agreement with several other models, such as the ones proposed by Laws 
et al. [23,24], and Nguyen [25]. 

From equation (7) it can be concluded that the plane stress compliance ten- 
sor of the damaged sublaminate 1, H (1 \ only contains two components that 
depend on the density of transverse matrix cracks: (L) and Hqq (L) . The 

tensor H (1 q is established as a function of the spacing of transverse matrix 
cracks, L, as: 


where: 


H (1 ) = 


' 1 
E~i 

V12 

Ex 

0 


V21 

e 2 

hS (i) 


0 

0 



h £ (i) = -=■ 


a£(l) 


d'’(i) A[\\L)Ag>(L) - \A f S(L) 


rd) 


<i), 


L<>), 


( 8 ) 


( 9 ) 


(i) 


1 

Gg(i) 


1 

^66'* {E) 


( 10 ) 


The evolution of the laminate stiffness as a function of the density crack 
f3 = 1/(2L) is shown in Figure 2 for a [±25/903] s laminate with the material 
properties shown in Table 1. The relation between the homogenized elastic 
properties of the cracked ply (sublaminate 1) and (3 is shown in Figure 3. 


Table 1 

Elastic properties of typical carbon/epoxy composite. 


E x (GPa) 

E 2 (GPa) 

G 12 (GPa) v \2 

163.4 

11.9 

6.2 0.30 
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Fig. 2. Laminate stiffness as a function of the matrix crack density in the 90° ply. 



Fig. 3. Effective elastic properties of the 90° ply as a function of the matrix crack 
density. 


Figures 2 and 3 show that the Young modulus of the laminate, E x , and the 
Young modulus of sublaminate 1, E \ , are not affected by the presence of 
transverse matrix cracks. Furthermore, the curves E 2 ((3)/E 2 and v 2 i(/3)/v 2 i 
are coincident (Figure 3), as a result of equation (7). 
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2.2 Onset of matrix transverse cracks under multiaxial loading 


The onset of transverse matrix cracking in a ply under the combined effect 
of in-plane shear stresses and transverse tensile stresses needs to be predicted 
using an appropriate failure criterion. The failure criterion should be estab- 
lished in terms of the actual strengths of a ply when it is embedded in a 
multidirectional laminate (in-situ strengths). 


2.2.1 In-situ strengths 

A methodology to predict the onset of matrix transverse cracking must be 
able to predict the in-situ strengths of a ply. Both the transverse tensile and 
in-plane shear strengths of a ply embedded in a multidirectional laminate 
are higher that the corresponding values measured in unidirectional laminates 

[31]. 

Models based on clastic analysis of the cracked laminate shown in Figure 1 
have been proposed to predict the in-situ strengths. Such an approach is only 
appropriate for the strength prediction of thin embedded plies. For a thin 
ply, it can be assumed that the defects at the micromechanical level can be 
represented as cracks that have a length equal to the ply thickness and can 
only propagate in the direction of the fibers [32], 

However, for a thick ply the length of the ’effective crack’ is smaller than the 
ply thickness and propagates in a first phase along the thickness direction [32] . 
Therefore, a methodology to predict the onset of matrix transverse cracking 
based on the cracked ply model represented in Figure 1 is not appropriate for 
thick plies. 

Instead, the thick ply model described in [32] is used. The tensile transverse 
in-situ strengths of sublaminate (1) is [32]: 


. / 8 Gi c 

For a thin embedded ply: i't = \ — — 

V 7 HA 2 2 

(ii) 

For a thick ply: Yf = 1.12\/2Yf d 

(12) 


where K^ d is the tensile transverse strength measured in a unidirectional test 
specimen, t is the thickness of sublaminate 1, G\ c is the mode I intralaminar 
fracture toughness, and A 22 is defined as: 



A° 2 = 2 


v. 


21 


E 2 E\ 


The in-situ shear strengths are obtained as [31]: 


(13) 


Sl = 


N 


(1 + x & G ^) 1 ^ ~ 1 

3x^12 


(14) 


where \ i s the shear response factor defined in [31], and the parameter <f> is 
calculated according to the configuration of the ply: 


For a thick ply: 


For a thin ply: 



(15) 

(16) 


where Sf d is the shear strength measured using an unidirectional test specimen 
and Guc is the mode 11 fracture toughness. 


2.2.2 Failure criterion for the prediction of transverse cracking under multi- 
axial loading 

In general, a ply represented by sublaminate 1 in Figure 1 is subjected to 
transverse tensile stresses and in-plane shear stresses. Under pure in-plane 
shear or pure transverse tension, the onset of transverse matrix cracking is 
predicted by comparing the components of the stress tensor with the respective 
in-situ strengths (defined in the previous section). 

Under multiaxial loading, it is necessary to use a failure criterion to predict 
the onset of matrix cracking. The LaRC04 [4,5] failure criteria are a function 
of the components of the stress tensor and in-situ strengths. For transverse 
tension, the criterion used is: 

+ + “ 1 - 0 with a22 - 0 ( 17 ) 

where g = Gi c /Gu c . Gn c is the mode II component of the fracture toughness 
associated with matrix transverse cracking. 
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Under moderate values of transverse compression, the crack plane is parallel 
to the laminate thickness direction. The LaRC04 failure criterion is: 

-r ( 1 0i2 1 + 7? L cr 22 ) - 1 < 0 with <r 22 < 0 (18) 

where r) L is the coefficient of longitudinal influence defined in [4,5]. 

2.3 Evolution of matrix transverse cracks under multiaxial loading 


To predict laminate failure under general loading, it is necessary to know how 
the density of transverse matrix cracks evolves in a ply subjected to combined 
in-plane shear and transverse tensile stresses. The evolution law is derived 
from the micromechanical model presented in this section. 

To define the damage evolution law, it is necessary to relate the applied stress 
or strain state to the density of transverse matrix cracks. This relation is 
obtained from a Fracture Mechanics analysis of cracked plies combined with 
the definition of the damaged constitutive tensor. 


2.3.1 Transverse tension 

Tan and Nuismer [18] obtained a closed-form expression that defines the evolu- 
tion of transverse matrix cracks under uniaxial stress states (either transverse 
tension or in-plane shear). 

To predict failure under multiaxial loading, i.e. when the lamina is simultane- 
ously under tensile and in-plane shear strains, £ xx and 7 xy respectively, it is 
necessary to derive a relation between the density of transverse matrix cracks 
and the applied multiaxial strain state. It is assumed that the relation between 
the tensile and shear strains is constant throughout the loading history: 

k = with e xx > 0 (19) 

£ XX 

where k is the multiaxial strain ratio. 

Consider a ply with crack spacing L, as shown in Figure 4 a). 
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a) Crack spacing 2L b) Crack spacing L 

Fig. 4. Progression of transverse matrix cracks. 

The strain energy in a laminate cell of length 2 L subjected to transverse 
tension and in-plane shear just prior to fracture, U 2 l, can be established as a 
function of the strain tensor and of the crack spacing as: 

U 2L — 2hLb E x (L)e 2 xx + A m (L)~{ 2 y (20) 

where E X (L) is the axial stiffness of the laminate with a crack spacing of 2 L 
and b is the specimen width; Aqq(L) is the laminate effective shear stiffness cor- 
responding to a crack spacing of 2 L, 2 h and 2 L are the laminate thickness and 
the distance between two consecutive transverse matrix cracks, respectively. 

After the propagation of transverse matrix cracks, Figure 4 b), the strain 
energy in the original unit cell of length 2 L is: 

U L = 2 hLb E x ^ e 2 xx + Aqq ^ rf xy (21) 

where E x (L/2) and Aqq(L/2) are respectively the axial stiffness and the lam- 
inate effective shear stiffness corresponding to a crack density defined by L. 

The energy required to generate a new matrix crack in a ply equals the loss of 
strain energy of the laminate [18]. Therefore the difference between equation 
(20) and equation (21) is equal to the energy released by the sublaminate 1: 

A U = U 2L -U L = 

= 2 hLb { - E, (1)] <4 + [/ML) - ^ (§)] 7^} (22) 


AU = 2h w bG c = 2h {l) bG ! + 2h {l) bG n (23) 
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where G c is the mixed-mode fracture toughness of sublaminate 1 under tensile 
(mode I) and shear (mode II) loading. From (22) and (23): 


E X {L ) — E x 


£ 2 + 


Am{L) ~~ A, 


66 


h^G r 


7 xy , 


hL 


(24) 


Using the definition of k given in (19), equation (24) can be re-written as a 
function of the strains: 


or: 


E X (L) — E x ( — ) ) + 


^66 (L) — A 66 ^ 


2 = h^G c 
hL 


(25) 


K- 


E X {L) — E x ^ + A m (L) — ^ | 


h^G r 


Ixy ~ 


hL 


(26) 


The relations between normal and shear strains, and the crack density are 
obtained as: 


&XX 


hWGr 


\ hL [E x (L)-E x ^)]+K 2 [A m (L)-A m ^) 


(27) 


7 xy 


hWG r 


K 


\ hL [e x (L)-E x ^)\+k*[A 66 (L)-A 66 (±) 


(28) 


Equations (27) and (28), are established as a function of the mixed-mode 
fracture toughness G c that needs to be defined. 

The mixed-mode fracture toughness is defined in terms of the mode I and 
mode II components as: 


G c = G\ + G a 


(29) 


From (29) and (24): 
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(30) 

(31) 


h^Gc 

hL 


h^Gi h^G u 
hL + hL “ 
E,(L) - E, (1) 


e 2 + 

^ XX ' 


Am(L) — A 66 




and: 


h^Gi 

hL 

h^Gn 

hL 


E,(L) - E„ (1) 
Am(E) — Aq$ (— ) 


£ 


2 

xx 



Dividing equations (32) and (33) by (25) and (26) respectively: 


(32) 

(33) 


A 




E X {L) — E x | 

( L \ 
\2j 


E X (L) — E x 

(f) 

+ K 2 

A 

66 (L) ~ A 66 | 

( L\ 
{ 2 ) 


A n = 


G\\ 

~Gc 


K 2 

Aee(L) — A 66 I 

f L \ 

l 


E X (L) — E x i 

(f) 

i 

+ K, 2 

Am 

(L) 

_ ^66 1 

( 1 ) 

1 


(34) 


(35) 


Ai and Ajj are the ratios between the mode I and mode II components of 
the energy release rate and the mixed-mode fracture toughness; Aj and An 
can be obtained as a function of the crack density and of the multiaxial strain 
ratio k. 


Taking into account that the energy release rate under mixed-mode loading is 
the sum of the mode I and mode II energy release rates, equation (29) can be 
established as: 


1 _ 

" G c + G c 

1 = At + A n 


(36) 

(37) 


The criterion proposed by Hahn [33] for the prediction of transverse matrix 
cracking under transverse tensile and in-plane shear loads is used: 


n 3 G He , , Gii , 

(! - 9) \ 7 ^ + 9 jr + 77 - = 1 
V Cric Oric LtHc 


( 38 ) 


Substituting equations (34) and (35) into (38) gives: 
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(1 ~9) 


1 


(39) 


IAjG c AjGc 
-Gi7 + 9 ^T 


+ 


Ic 


A n G c 

Guc 


or: 


(1 ~9) 


IA:G C 

G\c 


+ 9 


AjG c , (1 -Aj)G c 


G 


+ 


Ic 


G 


= 1 


He 


The positive real solution of (40) is: 


(40) 


G c 


G\\ c + 


41/ (Gic - Guc) 2 
2 Gle 



(41) 


where Aj, which is given by equations (34), depends on the density of trans- 
verse matrix cracks, /3, and on the multiaxial strain ratio, k (equation (19)). 

Figures 5 and 6 show respectively the relation between the crack density (3 
and e xx and 7 xy for different multiaxial strain ratios. A [±25/903] s laminate 
with the elastic properties shown in Table 1 is used in the predictions. 


£ xx ( mm / mm ) 



Fig. 5. Relation between applied strain (e xx and k) and density of matrix cracks. 
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Fig. 6. Relation between applied strain (7 xy and k) and density of matrix cracks. 

The effects of multiaxial strain states on the crack density are clearly shown in 
Figures 5 and 6. It can be observed in Figure 5 that the density of transverse 
matrix cracks increases with the multiaxial strain ratio for a fixed value of e xx . 


2.3.2 Transverse compression 

Transverse matrix cracks created by a combination of transverse compression 
and in-plane shear close under the effect of the compressive transverse stress. 
When a crack closes, its faces can transmit normal tractions but shear tractions 
may cause slippage between the crack faces. Therefore, it can be assumed that 
transverse matrix cracks only affect the shear stiffness of a laminate. 

Following the procedure described in the previous section, the energy released 
by the sublaminate 1 is: 


AU 


2 hL}ril y 


Am{L) — A 6e 



2h {1) bG Uc 


(42) 


The relation between the shear strain and the crack density is obtained as: 
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(43) 


7 xy 



hWG llc 

hi 

7^66 (L) — A 66 (7 


3 Micromechanics-based damage model 


3.1 Constitutive model 


A damage model able to represent the onset and accumulation of a periodic 
distribution of transverse matrix cracks should yield a compliance tensor sim- 
ilar to the one obtained from the micromechanical model, equation (8). An 
appropriate damage model can be developed by defining the Gibbs free energy 
per unit volume, Tg, as: 


1 

T c — -a : H : cr+ATa : cr + AM (5 : cr 
or, in the expanded form, as: 


(44) 


^G = 


a 


li 


+ 


a: 


22 


+ 


cr 


12 


E\ (1 — d 2 ) E‘2 (1 — do) G 


12 


^21 Vl2 

E'2 E\ 


V11V22 


+ 


+ (dllCTu + a 22 CT 22 ) AT + (/All'll + (^ 22 ^ 22 )AM 


(45) 


where Ei, E 2 , V\ 2 , v 2 \ and G\ 2 are the in-plane clastic orthotropic properties 
of a unidirectional lamina. The subscript 1 denotes the longitudinal (fiber) 
direction and 2 denotes the transverse (matrix) direction. d 2 and d$ are dam- 
age variables associated with transverse matrix cracking, an and 0:22 are the 
coefficients of thermal expansion in the longitudinal and transverse directions, 
respectively. and /3 22 are the coefficients of hygroscopic expansion in the 
longitudinal and transverse directions, respectively. AT and AM are respec- 
tively the changes in temperature and moisture content from the stress-free 
state. 

The proposed model distinguishes between active (d 2 +) and passive damage 
(g? 2— ) variables, corresponding to the opening or closure of transverse matrix 
cracks under load reversal respectively. The determination of the active dam- 
age variable can be accounted by the following equation: 
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( 46 ) 


d2 = d2+ <M +d2 _tM 
\^22\ \ a ’22\ 

where the McCauley operator (x) is defined as ( x ) := 2 (x + |x|). 

The damage variables d 2 and d 6 can be related to the density of transverse 
matrix cracks using the micromechanical model previously described. 

The constitutive model is obtained from the derivative of the Gibbs free energy 
with respect to the stress tensor: 

d^r 

e=— - = U : a+ATa + AM(3 (47) 

da 

where the compliance tensor H is defined as: 

0 

0 (48) 

1 

(1 — da ) G 12 _ 

The compliance tensor is established in terms of the damage variables is similar 
to the compliance tensor derived in the micromechanical model. It is important 
to note that the damage variables are applied to the H 2 2 and H m components 
of the compliance matrix because the micromechanical model results indicates 
that V2i{0) / E 2 {pi) = V 21 /E 2 is independent of damage (Figure 3). 

Using the stiffness tensor C = HU 1 , it is possible to re-write equation (47) as: 

a = C : (e-A Ta - AM/?) (49) 

with: 

Ei (1 — d 2 ) E 2 V 12 0 

C —— (1 — d 2 ) E 1 V 21 (1 — d/ 2 ) E 2 0 (so) 

0 0 A (1 - d 6 ) G 12 

and A = 1 - u 2 Uh 2 (l - d 2 ). 


V21 

Ei E2 

d 2 T g _ V12 1 

da 2 ~ Ei (1 - d 2 ) E 2 

0 0 
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In the context of Continuum Damage Mechanics it is necessary to relate the 
effective stress tensor, a, used in the damage activation and evolution functions 
to the the nominal stress tensor a. The relation between the effective and 
nominal stress tensors is established in the proposed model using the Principle 
of Strain Equivalence [34]: 


a = M : a (51) 

where the damage operator, M, can be obtained from equations (48) and (50): 


M = CU =0 : H 


1 ^ 12^2 

(1 - d 2 ) (1 - V 21 V 12 ) 

q 1 — ^ 21^12 + V 2 iZ>i 2 d 2 

(1 - d 2 ) (1 - v 2 iv 12 ) 
0 0 


0 

0 

1 

1 — d,Q , 


(52) 


Using (8) and (48) the damage variables d 2 and d$ can be expressed in terms 
of the crack density (3 as: 


d'2+ — 1 — 


1 

E 2 H 22 {(3) 


(53) 


dg — 1 — 


1 

G\ 2 Hqq (/ 3 ) 


(54) 


The functions H 22 (/ 3 ) and H 66 (/ 3 ) are obtained from equations (9) and (10) 
with (3 = 1/(2 L). 

It can be shown that the damage variables are defined between zero, when 
crack spacing tends to infinity, and one, when the crack spacing tends to zero. 
The relation between the density of transverse matrix cracks f3 and the damage 
variables d 2 and d 6 is shown in Figure 7. 


18 




Fig. 7. Relation between the damage variables d 2 and d§ and density of transverse 
matrix cracks. 

The condition of positive energy dissipation is established as a function of 
the thermodynamic forces, Y, and as a function of the time derivative of the 
damage variables, d, and it is given by: 


or 


S = Y d 


d^ G 

dd 


■ d > 0 


o. 


22 


2(1 — G?2+) E 2 


do+ + 


cr 


12 


2 (1 — cIq) Ct 12 


~d( 3^0 


(55) 


where S is the rate of energy dissipation per unit volume. 

The rate of dissipation can be established in terms of micromechanical vari- 
ables by using equations (46-54): 


1 

(1 — C ? 2+) 2 E 2 


[H 22 (I3)} 2 E 2 


(56) 


1 

(1 — do) 2 G12 


\H m m 2 o 12 


The time derivatives of the damage variables are: 


( 57 ) 
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1 


(58) 


d'2+ — 


d'6 


1 dH 22 (P) 

\H 22 {(3)fE 2 dp 

1 dH m {P) 


P 


[H m ((3)} 2 G 12 dp 
Therefore, the rate of energy dissipated is: 


P 


(59) 


9H 22 (P) ^ + a\ 2 dH m (P) 


dp 


dp 


(60) 


Taking into account that dH ^^ > 0 and dH g 6 S ^ >0, the time derivative of [3 
must be greater than or equal to zero in order to assure positive dissipation, 
he., P> 0 . Therefore, the condition of positive dissipation, when interpreted 
from a micromechanical point of view, establishes that the density of trans- 
verse matrix cracks can only increase or remain constant. 


3.1.1 Damage activation functions 

Transverse matrix cracks are predicted using two scalar functions, F k (k = 
2+, 2—), established in terms of the effective stress tensor a*, and of the dam- 
age threshold value, r l : 


Fk <j>k (cr*) - r* < 0 (61) 

where t is the current time. The initial threshold value, r°, is equal to 1, and 
the following condition must be satisfied in order to fulfill the Second Principle 
of Thermodynamics: 


f > 0 (62) 

Damage onset occurs when any of the functions (p k (d 4 ) reaches the initial 
damage threshold value of 1. The functions (p k (d f ) used to predict transverse 
matrix cracking are based on the LaRC04 failure criteria previously proposed 
by the authors [4,5]. 

The damage activation functions are used to predict the onset of transverse 
matrix cracks lying in ply thickness direction, as shown in Figure 1. 
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Transverse matrix cracks lying in the direction of the ply thickness occur 
under transverse tension (<r 2 2 > 0), or under moderate values of transverse 
compression and in-plane shear (cr 22 < 0). 

For high values of transverse compression, the matrix fractures lie along a 
plane that is inclined at an angle a to the direction of the ply thickness, as 
shown in Figure 8. 



Fig. 8. Matrix crack in a created by high in- plane compressive transverse stress. 

Inclined fracture planes caused by transverse compression induce delamination 
between the plies. This damage mechanism is usually catastrophic in uniformly 
stressed composites where local redistribution to more lightly loaded regions 
of the structure cannot occur [6]. Therefore, laminate failure can be assumed 
to occur when matrix cracking under high values of transverse compression is 
predicted. 

The damage activation function used to predict matrix cracking under trans- 
verse tension (at 22 > 0) and in-plane shear is defined as: 


F 2+ := 0 2+ (d 4 ) — r 4 < 0 


(63) 


with: 


02+ (cf 4 ) 


M 


a - 9 ) 


y t 


+ 9 




(64) 


The damage activation function used to predict matrix cracking with a = 0° 
under moderate values of transverse compression (d 22 < 0) and in-plane shear 
is defined as: 


F 2 - := 02- — r 4 < 0 


(65) 
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with: 



cr 


12 


+ r ] L a l 


22 


( 66 ) 


The equations proposed can accurately predict transverse matrix transverse 
cracking. The failure envelope, or initial elastic limit, predicted using equations 
(63) and (65), and the corresponding experimental data are shown in Figure 
9. 



Fig. 9. Comparison between predictions and experimental data [36] and definition 
of ply or laminate failure domains. 

Figure 9 also shows the domain of validity of the model: the damage model is 
only defined for fracture angles a = 0°. For a > 0° the laminate is assumed 
to fail. 


22 


3.1.2 Damage evolution functions 

Using the principle of strain equivalence, the effective stress tensor can be 
written as a function of the strain and undamaged stiffness tensors: 


a = C°: e (67) 

The strain tensor is established in terms of the density of transverse cracks 
and the multiaxial strain ratio k, equations (27-28). Therefore, all terms in 
the damage threshold functions can be formulated as a function of the density 
of transverse matrix cracks and of the the multiaxial strain ratio: 


F t = f 0?V) 


( 68 ) 


For a given state of strain at time t the multiaxial strain ratio is a dependent 
variable that can be easily defined using equation (19) in material coordinates: 




(69) 


The density of transverse matrix cracks is a state variable. Therefore, it is 
necessary to define an evolution law subjected to thermodynamic restrictions. 

The first condition to be satisfied is the requirement of positive dissipation. As 
demonstrated by equations (56)-(60), the condition of positive dissipation is 
satisfied when the evolution of the state variable f3 is defined by a monotonic 
increasing function: 


/3 > 0 (70) 

Furthermore, it is necessary to define the conditions for the evolution of the 
clastic domain: 


Fkia*,/?) <0 

(71) 

$F k (a t ,0 t ) = 0 

(72) 


Equations (70)-(72) are the Kuhn- Tucker conditions which ensure a consistent 
evolution of damage during loading and load reversals. 

The evolution laws of the state variables are defined as: 
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(73) 


da . da 
a de de 


dr ■ dr 

r= aU + ^ K 




C° : 




(74) 


For a given loading state, the damage consistency condition must be applied 
to define the evolution of the internal variables. The consistency condition is 
defined as: 


dF k • dF k 

F k = 0 => F k = — - : a H — —^r = 0 
oa or 


(75) 


From (63) and (65): 


dFk 

dr 


= -1 


(76) 


Using (76) in (75): 


• dF k ■ . dF k . . 

F), = - : a — r = - 1 ^ r : C : £ — r = 0 




<9d 


(77) 


Using equations (73) and (74) in (77): 


6-S^(S»*l-)-(S **!*)-» « 


Taking into account that the micromechanical model proposed assumes a con- 
stant loading ratio, k = 0: 


F k 


dFk ( C o.fe 

da ' { ' dp 



(79) 


The density of of transverse matrix cracks, /?, is calculated from the integra- 
tion of equation (79) using a numerical method, such as the return-mapping 
algorithm. 

For damage evolution under transverse tension, the values of dF k /da can be 
obtained from (63): 


dF 2+ 

dan 


0 


(80) 
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1 l-g 1 


^22 


( 81 ) 


dF 2+ 
da 2 2 


02+ 


+ Q 

(Ft) 


dF 2+ = 1 CTl2 

C>ai2 02+ (S L ) 2 


and for damage evolution under transverse compression (65): 


(82) 


dF 2 . 
da xl 


= 0 


(83) 


dF 2 . 

da 22 


V_ 

Sl 


(84) 


dF 2 . _ I ^-di 2 if a 12 > 0 
d<Tl2 \ if a V2 < 0 


St 


Cl2 


(85) 


Figure 10 represents the evolution of the elastic domain for different values of 
the crack density f3 using a ply with the material properties shown in Table 1. 



Fig. 10. Predicted evolution of the elastic domain as a function of crack density (3. 
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4 Prediction of laminate failure 


In the presence of stress concentrations the onset of fiber localized failure 
does not cause immediate structural collapse. Experimental results [37,38] 
have shown that structural collapse is caused by the the progression of fiber 
fracture. Therefore, it is necessary to use a damage model that accounts for 
the stress re-distributions caused by fiber fractures. 

The model proposed here assumes that for uniformly stressed laminates the 
onset of fiber fracture and delamination caused by high compressive trans- 
verse stresses triggers structural collapse. Therefore, laminate final failure is 
predicted when fiber failure or matrix cracking with « ^ 0° occurs. 


f.l Fiber failure 


The criterion for fiber fracture under longitudinal tension (<7n > 0) is defined 
as [4,5]: 


FIi+ := A — 1 < 0 (86) 

Ax 

where is the ply tensile strength in the longitudinal direction. 

The LaRC03 failure criterion for fiber kinking under longitudinal compression 
(<r ii < 0) is a function of the components of the stress tensor in a frame 
representing the fiber misalignment, a [4,5]: 


(Jif* = an cos 2 ip + <7-22 sin 2 (f + 2 \u 12 1 sin ip cos ip 

ofn = (Tn sin 2 ip + a 22 cos 2 ip — 2 |<ti 2 | sin ip cos ip (87) 

<j 12^ = —an sin ip cos p> + (T 22 sin ip cos ip + |(Ti2 1 (cos 2 ip — sin 2 ip^j 

where the fiber misalignment angle ip is defined as [4]: 


<P 


012 1 + (Gl2 ~ Xq ) ip C 

G 12 + (Til — <^22 


( 88 ) 
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where 


ip c = tan 1 



(89) 


and where Xc is the ply compressive strength in the longitudinal direction. 

Depending on the sign of the in- plane transverse stress (J%f\ the criteria for 
fiber kinking (<7n < 0) are: 


FIi_ : = 


M m ) 

°12 


+ yd?’ 


S,. 


- 1 < o, yy < o 


(90) 


or 


FIi- : = 


(m) / (m ) ' 

\ a i2 i | , 

(! -9) ~^+9 I ww I + 


a 12 


M' 


Kt 


It 




- 1 < o, 4’ 2 n) > 0 


(91) 


4-2 Matrix failure with a ^ 0° 


The failure criteria for matrix cracking under transverse compression (<7 2 2 < 0) 
and in-plane shear and a ^ 0° are defined as [4,5]: 


Fla. := 




- 1 < O,0ii > — Y c 


(92) 


FI 2 _ := 



2 

- 1 < 0,0-n < -Y c 


(93) 


where the effective shear stresses in the fracture plane are defined as: 


tT — 


t l = 

' a 


T 


+ T] T <7 n COS O^j 

+ rj L <7 n sin O'j 


(94) 


where 9 = tan 1 ( — \^L ) . The components of the stress tensor on the fracture 

V 0-22 sin a J L 

plane are given by [4,5]: 


27 



< 


( 95 ) 


o n = 0 2 2 cos 2 a 
t t = — o 2 2 sin a cos a 
t l = o 12 cos a 

The terms rJ" T and t ( T l are calculated from equations (94)- (95) using the 
relevant components of the stress tensor established in a frame representing the 
fiber misalignment, equation (87). The angle a is determined by maximizing 
the failure index FI 2 - (92-93) using a simple iterative procedure. 

The coefficients of transverse and longitudinal influence, 77 T and rj L respec- 
tively, are [4,5]: 


T 1 

77 = 

tan 2a 0 

(96) 

L Sl cos 2a 0 

n = v 2 

Y c cos- ao 

(97) 


with ao ~ 53°. Yq is the ply compressive strength in the transverse direction. 


5 Examples 


The present damage model can be used in combination with classical lamina- 
tion theory using stand-alone codes. Alternatively, the damage model can be 
implemented as a constitutive subroutine for the Finite Element Method. 

The damage model was implemented using a commercial symbolic computing 
software. The model was verified by calculating the response of several glass- 
epoxy laminates under uniaxial and multiaxial loads: [±45/904] s , [02/904] s , 
[0 2 /90 2 ] s and [0 2 /90] s . 

In all calculations performed, the ply thickness was taken as 0.144mm and the 
temperature difference from the stress free condition -100°C. The coefficients 
of thermal expansion in the longitudinal and transverse directions are an = 
7.43 x 10 _6 /°C and 022 = 22.4 x 10 _6 /°C, respectively. The remaining material 
properties used are shown in Tables 2 and 3. 
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Table 2 


Elastic properties of glass-epoxy [20] 


Ei (GPa) 

E 2 (GPa) 

G12 (GPa) 

G 23 (GPa) 

V 12 V 23 X (10 8 MPa 3 ) 

44.7 

12.8 

5.8 

4.5 

0.30 0.42 2.0 


Table 3 


Strengths and fracture toughnesses of glass-epoxy [20] . 


Y| d (MPa) 

Sl d (MPa) 

Gi c (N / mm 1 ) 

Guc (N/rnrn x ) 


40.0 

73.0 

0.20 

0.40 



The predicted response of [±40/904] s and [02/904] s , is shown in Figure 11 for 
two values of the multiaxial strain ratio: k = 0 and k = 10. 

Ex(p)/Ex 



s xx (mm /mm) 


Fig. 11. Relation between laminate modulus and applied strain for [±45/904],, and 
[O2/9O4] s laminates. 

Figure 11 shows that the rate of degradation of the elastic properties of the 
laminate is higher when the axial stiffness of the outer sublaminate decreases. 
The effect of multiaxial loading is also clear in Figure 11: as expected, the 
application of shear strains leads to a reduction of the extension corresponding 
to the onset of transverse matrix cracks and to a higher rate of degradation 
of the elastic properties of the laminate. 
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Figure 12 compares the response of [ 02 / 904 ) 5 , [O 2 / 9 O 2 ] s and [O 2 / 90] s laminates 
for k = 0 and k — 2. 


Ex(|3)/Ex 



e xx (mm/mm) 


Fig. 12. Relation between laminate modulus and applied strain for [ 02 / 904 ) 5 , 
[O 2 / 9 O 2 ] s and [02/90] s laminates. 

The in-situ effect is shown in Figure 12. For k — 0, the strain corresponding 
to the onset of matrix cracking of the [02/90] s laminate is 1.9 and 1.3 times 
higher than the ones of the [0 2 /90 4 ]s and [0 2 /90 2 ]5 laminates respectively. 
Furthermore, the strain at the onset of matrix cracking of the [O 2 / 90] s is 2.7 
times higher than the ultimate transverse strain measured in an unidirectional 
test specimen. 

The laminate modulus calculated using the proposed model and the ply dis- 
count method is shown in Figure 13 for a [0 2 /90 4 ]5 laminate. The ply discount 
method assumes that E 2 ~ 0 and G 12 ~ 0 as soon as the failure criterion is 
satisfied. 
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Ex(p)/Ex 


(0,/90. 


2 ' i's 


K— 0 


1 , 0 - 


0,9- 


0 , 8 - 


Damage model 


0,7- 


Ply discount method 


0,6 -| 1 1 1 1 1 i i 1 i 1 1 1 1 1 

0,000 0,002 0,004 0,006 0,008 0,010 0,012 0,014 

e xx (mm/mm) 


Fig. 13. Laminate modulus calculated using the damage model and the ply-discount 
method. 

It can be observed in Figure 13 that the response of the modulus of the 
laminate obtained using the ply-discount method is a step-function that does 
not accurately represent the progressive degradation of the laminate. 


6 Conclusions 


A new, micromechanics-based, continuum damage model able to simulate the 
onset and propagation of transverse matrix cracks and final laminate failure 
is proposed. The model is applicable to [±0/9O n ] s laminates, under multiaxial 
loading and uniform stresses or small stress gradients. 

The model uses ply properties and does not require any tests performed at 
laminate level to identify damage onset and evolution functions. The onset 
of damage is predicted using failure criteria and damage evolution laws are 
established from the micromechanical analysis of cracked plies. 
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The model can be used in Finite Element analysis, and the results are inde- 
pendent of the mesh refinement because the constitutive model does not result 
in strain-softening. The onset and accumulation of transverse matrix cracks 
are represented as a distributed damage mechanism. The onset of localization, 
which is triggered by either fiber fracture or matrix cracking with a ^ 0, is 
assumed to cause a structural collapse. 

The predictions show that the rate of degradation of the clastic properties 
increases when the stiffness of the outer sublaminate decreases. Decreasing 
the thickness of the 90° plies also increases the degradation rate. 
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Appendix A: Laminate constitutive model 


The plane stress constitutive equation of a laminate is: 

£ = Q _1 a (A-l) 

where the barred parameters denote a quantity that is not a function of the 
thickness coordinate. 

The laminate stiffness Q matrix can be obtained from the stiffness and thick- 
ness of the individual plies applying the classical in-plane laminate theory: 

i n 

Q = rX> (fc) Q (fc) (A-2) 

' l k=\ 

where h is the total laminate thickness, hS k ' 1 are the thickness of the k th ply and 
Q is the stiffness of the each k tn ply in global laminate coordinate system 
(see Figure 1) 

For a balanced symmetrical laminate the matrix Q is defined as: 

Qn Q12 0 

Q = Q21 Q22 0 (A-3) 

0 0 Qqq 

The stiffness of a laminate composed by two sublaminates is obtained as: 

Q = -L (2/r (1 )Q« + 2h^Q {2) ) = i (/r (1) Q (1) + /i (2) Q (2) ) (A-4) 

For a balanced, symmetric laminate, the laminate stiffness depends on the 
distance between two consecutive longitudinal cracks (L): 

A(L) n A(L) 12 0 

A (L)= A(L) 21 A(L) 22 0 (A-5) 

0 0 A (/) 66 

where A lJ (L) 1 i, j = 1,2 are the average laminate stiffness: 
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A n {L) 


+ h^Q ( S 

(hh 


(A-6) 


Ai2(L) 


P 2 h 


(A-7) 


A22 (A) 


/i (1) Q& + hWQ$ 
h 


( fc-pA^ fcg) 2 

V & ) hQn 


(A-8) 


A m (L) is the laminate shear stiffness calculated as: 


A-mP) 


h'VQtt + hWQ® 

@4 h 


(A-9) 


The superscript (1) and (2) denote ply 1 and ply 2 respectively and = 

1, 2, 6 are the undamaged lamina stiffness along the laminate axes. 

The parameters 0i and 0 2 can be calculated as [18]: 


A = 1 


02 = 1 + 


tanh a \ L 
otiL 

(A-10) 

h^Qu \ tanh a\L 

h ( 2 '>Q ( u> I ay L 

(A-ll) 


with ay: 


2 = 3 Cgcg ( h^Q[\ } +h^Q ( pA 

hWC$ + hWcjg V hmWQftQ?! ) 

where CiP is the lamina stiffness in the lamina axes. 

The parameter 0 4 is be calculated as: 


h { 1 ) Qm\ tanh a 2 L 

h'' ,, Q'S) L 


(A- 13) 


with a 2 : 
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(« 2 ) 2 


3 C$c£> / fe (1) g^ + fe (2) gg } \ 

hWC$ + hWC$ \ hWhWQ$Q$ ) 


(A- 14) 


where C]$ is the lamina stiffness in the lamina axes. 

The relation between the longitudinal modulus E x , the Poisson ratio v xy and 
the shear modulus G xy and the crack density is obtained Aij(L),i,j = 1,2,6. 

After the definition of Aij(L),i,j = 1,2,6 it is possible to calculate the en- 
gineering constants of the laminate. The longitudinal modulus, E X (L), of the 
quarter cell under plane stress is: 


E x (L) = A n (L ) 


Au(L) 

A 22 (-£') 


The Poisson ration of the laminate is: 


v xy(L) 


AuiL) 
A22 (L) 


The shear stiffness of the laminate, G xy (L), is obtained as: 


(A-15) 


(A-16) 


G X y(L ) — Aqq(L) 


(A-17) 
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